Breakdown of self-averaging in the Bose glass 
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We study the square-lattice Bose-Hubbard model with bounded random on-site energies at zero 
temperature. Starting from a dual representation obtained from a strong-coupling expansion around 
the atomic limit, we employ a real-space block decimation scheme. This approach is non-perturbative 
in the disorder and enables us to study the renormalization-group flow of the induced random-mass 
distribution. In both insulating phases, the Mott-insulator and the Bose glass, the average mass 
diverges, signaling short range superfluid correlations. The relative variance of the mass distribution 
distinguishes the two phases, renormalizing to zero in the Mott insulator and diverging in the Bose 
glass. Negative mass values in the tail of the distribution indicate the presence of rare superfluid 
regions in the Bose glass. The breakdown of self-averaging is evidenced by the divergent relative 
variance, increasingly non-Gaussian distributions, and a correlation-length exponent v = 0.6 ± 0.05 
that violates the Harris-criterion bound. 
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Introducing quenched disorder into an otherwise pure 
system can lead to subtle and complex results and the 
disordered Bose-Hubbard (BH) model is no exception. 
While the clean BH model shows a relatively straightfor- 
ward bosonic competition between repulsion and tunnel- 
ing, the disordered model exhibits a new gapless insulat- 
ing phase, the Bose glass (BG), the precise location of 
which has proven problematic from the outset (T]. With 
the advent of experimental methods that can engineer 
this model directly [2J, the problem has been inverted 
and this has sparked renewed interest in the role disor- 
der plays in quantum systems. 

In this Letter, we develop a non-perturbative method 
to probe the nature of the transition between the 
localization-induced BG and the Mott insulator (MI) in 
the two-dimensional BH model with bounded potential 
disorder. The MI arises from the on-site repulsions and 
hence dominates in the limit the hopping vanishes while 
the superfluid (SF) is the ground state in the opposite 
regime. It is in the difficult intermediate parameter space 
where the BG phase obtains. Fisher and colleagues p] 
argued that the BG is a quantum Griffiths phase domi- 
nated by arbitrarily large SF regions that are, however, 
exponentially suppressed. Despite the abundance of nu- 
merical [3HS] and analytical [TCJH2CJ] work on the subject, 
it is only recently that several aspects of this model have 
been fully understood, such as confirmation that the BG 
always intervenes between MI and SF phases [2T] (see 
Fig. [IJ and the distinction between the MI and BG re- 
garding whether fluctuations are self-averaging [22] , 

As argued by Aharony and Harris [23], the break- 
down of self-averaging can be identified from the 
renormalization-group (RG) flow of the relative variance 
of any extensive variable. If the relative variance does not 
renormalize to zero, the central limit theorem no longer 
applies and the system is not self-averaging. This con- 
cept has been used to characterize the phase transition 




FIG. 1: Sketch of the phase diagram for the disordered Bose- 
Hubbard model as a function of the hopping t/U and chemical 
potential fi/U in units of the on-site repulsion U. Random 
shifts in the on-site energies bounded by A < U/2 lead to the 
formation of a compressible Bose-glass (BG) phase, separat- 
ing the Mott-insulating (MI) lobes from the superfluid (SF). 



between the MI and the BG [22] within a disorder av- 
eraged Replica field theory which is valid close to the 
atomic strong-coupling limit. In both insulating phases, 
the mass of the theory diverges, signaling the presence of 
short-ranged SF correlations. Since in dimensions d < 4 
the variance of the mass distribution diverges as well, the 
breakdown of self-averaging can be readily understood 
as a competition between the spread of the distribution 
versus the shift of its average. In the MI, the shift domi- 
nates the spread leading to a vanishing relative variance. 
In the BG, the spread is faster and the relative variance 
diverges. This characterization of the MI to BG tran- 
sition is consistent with the commonly accepted picture 
that the BG is dominated by rare SF regions since the 
negative mass values occur in the tail of the distribution. 
Whether the onset of this Griffith-type instability is cor- 
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rectly described by an RG calculation that is perturbative 
in the disorder remains a central question. 

We address this question within a real-space block dec- 
imation RG scheme on a d = 2 dimensional square lat- 
tice. This approach is non-perturbative in the disorder 
and hence enables a study of the RG flow of the full 
random-mass distribution, a necessity for any definitive 
statement about Griffiths-type physics to be made. Our 
results confirm that the relative variance serves as the or- 
der parameter for the MI-to-BG transition. Determining 
the correlation length from the scale at which the relative 
variance becomes of order one, we extract a correlation- 
length exponent v = 0.6 ± 0.05 which is close to the 
analytical value v = X/d obtained within the perturba- 
tive 1-loop RG [22 . It has been argued [23l EI] that a 
violation of the Harris-criterion bound v > 2/d for crit- 
ical disordered systems E5] EE] is indicative of the lack 
of self-averaging. The absence of the central-limit the- 
orem in the BG is further evidenced by an increasingly 
non-Gaussian shape of disorder distributions. 

Our starting point is the simplest form of the disor- 
dered BH model on a square lattice, 

H = -tJ2 (blbj + fttfei) - IHfk + -g £ *i (»i - 1) , 

(1) 

which describes bosons tunneling with amplitude t be- 
tween nearest-neighbor sites i and j and interacting via 
an on-site repulsion U. The bosonic raising and lowering 
operators are given by b\, bi respectively where hi = b\bi 
is the bosonic number operator. Hi is the chemical poten- 
tial shifted by the on-site disorder potential, Hi = h — e;. 
The random site energies are uncorrelated between 
different sites and uniformly distributed in the interval 
[—A, A]. From minimization of the energy in the atomic 
limit it is straightforward to see that for A < U /2 the 
phase diagram retains MI phases (see Fig. [T]) . 

To facilitate a strong coupling expansion we follow the 
standard procedure E7] EH] . After expressing the model 
by a coherent-state path integral in imaginary time, we 
decouple the hopping term by a Hubbard-Stratonovich 
transformation and trace over the original boson fields. 
We then perform a temporal gradient expansion to obtain 
the effective dual action on a lattice, 



S e ff = a' 



(2) 



where a denotes the lattice spacing and the sum over 8 
runs over the nearest neighbors of site i. By construction, 
the complex fields correspond to the SF order parameter 
ipi ~ (bi)- In regions where Ri > 0, SF order is sup- 
pressed. The mass Ri therefore corresponds to the local 
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FIG. 2: RG flow of the mass distribution P{ri/rl) normalized 
to a mean of 1 in the MI. Starting with the initial distribution 
given by Eq. (3a I, each subplot represents the probability 



distribution of this quantity after another iteration of the RG 
equations. We see that the Gaussian fit (red dotted line) 
becomes better with successive iterations. In addition, the 
width of the distribution, which corresponds to the relative 
variance, narrows. 



Mott gap. The calculation of the coupling constants as a 
function of the microscopic parameters of the BH model 
is straightforward, e.g. Ri and Hi are related to the 
single- and two-particle bosonic Green functions of the 
local on-site Hamiltonian, respectively, while the tempo- 
ral gradient terms and -JQ are given by derivatives 
of the mass with respect to the chemical potential. We 
specialize to the first Mott lobe with m — 1 bosons per 
site in which the coefficients are given by []j|J EH] 
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where z — 2d — 4 is the coordination number of the 
square lattice. In the clean limit, Hi — Mi f ne mean-field 
phase boundary between the first MI lobe and the SF is 
obtained from R(h, t, U) = 0. In the presence of disorder, 



Mi 



//—£,, the coefficients Ri, -K,- , and Hi depend 



on the disorder potential e,, which induces non-trivial dis- 
order distributions of the coefficients. Note that initially 
the dual hopping amplitudes T,<5 = l/{za 2 t) are uniform. 
We allow for a spatial dependence of the hopping since 
disorder will be induced under block decimation. 

Equations ^ provide the initial conditions for our pro- 
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cedure where the random on-site energies are generated 
from a uniform distribution on the interval [—A, A]. In 
order to determine the phase for a given set of parameters 
t, U, /i, and A of the disordered BH model ([T]), we derive 
a set of recursion equations using block decimation. 

We start with the discrete action ^ and eliminate 
short-range degrees of freedom by integrating out every 
other site, treating the quartic terms Hi perturbatively to 
leading order. The remaining points form a new square 
lattice with lattice spacing a' = \[2a and tilted 45° from 
the original system. A second iteration of this procedure 
brings the system back to the original angle with appro- 
priate reindexing. Recollecting the resulting terms and 
rescaling the action to look like the original, we find the 
RG recursion equations to be 



flj, = Ri + Tis - T iS T iS'h (l - misle) (4a) 

S 5,8' 

T h = J2 T eTe>Iee> (l - 4H ee , I ee , I e£ , ) , (4b) 

where Ii = (Ri + Tig) -1 is the static propagator and 



(4(iZi + E 



S T U )K 



(2) 



iX)\2 



The indices 5 



and 6' correspond to nearest neighbors of site i, whereas 
e and e' correspond to the bonds adjacent to the bond 
connecting the sites i' and f of the remaining lattice. 
The site ee' is the common vertex of the bonds e and e' 
(see Fig. [3]). Since we are interested in the MI to BG 
transition at incommensurate filling, we can neglect any 
corrections to the coefficients K^ 2 \ and H beyond 

dimensional scaling. The dynamic exponent is d z = 2 
since ^ except for a vanishing subset of sites of 
zero measure. 




FIG. 3: Illustration of the block decimation scheme underly- 
ing the recursion relations Q. 



From Eq. (4b) it is obvious that the gradient terms 
are renormalized under block decimation. The effective 
mass should therefore be defined relative to the kinetic 
hopping terms. In the following, we will use the effective 
local mass variable 
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FIG. 4: RG flow of the mass distribution P(ri/ri) in the BG. 
Note the similarity to Fig. [2] up to an overall displacement 
in the first histogram. Here we see the Gaussian fit becomes 
worse and worse due to the long tail of the distribution. This 
corresponds to an increasing relative variance. 



As a consistency check, we evaluate our recursion re- 
lations in the clean limit. For T — Tis, R = Ri, 
and H = Hi = 0, the RG equations Q reduce to 
R! = R + AT - 16T 2 /(R + AT) and T = 47*/(i? + 4T), 
leading to the recursion relation r' = 2r + r 2 /A for the 
effective mass. This indeed correctly describes the mean- 
field transition between the MI and the SF at r — 0. 

In the following, we integrate the RG equations Q 
numerically for the inhomogeneous system obtained for 
one particular disorder realization and keep track of the 
values of the coupling constants on each lattice site. 
This allows us to extract the mass distribution P(ri) at 
each iteration step of the block decimation. While we 
vary the hopping t/U and the lattice constant a to tune 
through the MI-to-BG transition, we use the same val- 
ues A/U = 0.1 for the disorder strength and [ij\J = 0.15 
for the chemical potential throughout the paper. Since 
Eq. ([2]) has only nearest neighbor or on-site terms, we can 
ignore the boundary points after each iteration without 
affecting the overall distribution. The major limitation 
to this method is the necessity of finite size lattices. Es- 
timates of any diverging quantities near the critical point 
must take into account finite system size effects. Accom- 
panying this limitation is the computation time increase 
for correspondingly larger system sizes. In the data given 
below we use an initial square grid of points with side 
length L — 506 sites. 

As expected for insulating phases, in both the MI and 
the BG, the mean r% of the mass distribution increases 
exponentially under the RG signaling short-range SF cor- 
relations. To distinguish the behavior in the MI and the 
BG, we normalize to a mean of unity and analyze the 
evolution of P(ri/rt). The variance of this rescaled dis- 
tribution corresponds to the relative variance of the mass 
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FIG. 5: Scale collapse of the correlation length plotted as a 
power law of the reduced hopping. The fitting parameters 
are t c and v. This fit corresponds to the least square error 
confirming v — 0.6 ± 0.05. 



distribution, which should serve as order parameter |22j . 

In Fig. [2] the RG flow of P(r l /f^) in the MI phase is 
shown. Note that the initial mass distribution is asym- 
metric and non-Gaussian due to the functional depen- 



dence (3a I on the uniformly distributed on-site energies 
Ci . As a consequence of the relatively small hopping value 
t/U, the bulk of the sites begin well above = and 
continue together towards larger values under repeated 
iteration. Increasing the number of iterations results in 
a distribution well described by a Gaussian. Further, the 
width of the distribution narrows, indicating a vanish- 
ing relative variance. This demonstrates that in the MI 
disorder is irrelevant and the system is self-averaging. 

The situation changes dramatically as we increase the 
value of t/U and enter the BG phase (see Fig. [4]) . While 
the overall shape of the initial distribution looks quite 
similar to the one in the MI, now a large fraction of the 
initial sites lies close to or below — 0. With enough 
sites close to the transition point, some regions of the sys- 
tem take much longer to flow to large positive values than 
others. This results in a drastic spread of values upon 
repeated iteration of the RG, leading to a divergence of 
the relative variance. In addition, the distribution de- 
velops a non-Gaussian form under successive iterations, 
indicating a violation of the central limit theorem. This 
demonstrates a breakdown of self-averaging in the BG. 

Our results show that the relative variance of the 
random-mass distribution serves as an order param- 
eter for the transition between the MI and the BG. We 
can therefore determine a correlation length £/a = \[2 
in the BG from the average number of iterations it takes 
before the relative variance becomes of order unity. For 
the parameters used in Fig.[4]this happens between 6 and 
7 iterations. A more precise value is obtained by aver- 
aging over several disorder realizations. Note that once 
the relative variance becomes of order unity, the left tail 



of the distribution pushes through zero. Therefore, the 
correlation length £ corresponds to the typical distance 
between SF droplets in the system. 

At the transition to the MI we expect the correlation 
to diverge as a power law £ cx (t — t c )^' / . To extract 
the correlation-length exponent i>, we vary the hopping 
slightly above the transition point for a fixed value of 
the initial lattice spacing a and extract the correlation 
length as described above by averaging over several dis- 
order realizations. In the following, we use lattice con- 
stants a = 0.2,0.25,0.3,0.35,0.4 and average over 10-20 
disorder realizations for each value. Note that the value 
t c of the transition point is non-universal and depends 
on a. As shown in Fig.[5j plotting £/a vs [t — t c {a)]~ y /a, 
the different data sets collapse on a single straight line 
for v = 0.6 ± 0.05. This value violates the Harris cri- 
terion bound v > 2/d [251 ES] but is consistent with a 
breakdown of self-averaging [23] HI] ■ 

Discussion and Conclusion. We used a real-space 
block decimation method to characterize the MI-to-BG 
transition of the disordered BH model. By analyzing the 
RG flow of the induced random-mass distribution in the 
dual field theory, we have demonstrated that the tran- 
sition is characterized by a breakdown of self-averaging. 
The associated correlation length corresponds to the typ- 
ical separation of rare SF regions. Our work provides the 
first explicit confirmation that the instability towards the 
formation of the BG is of the Griffith type. 

The method employed here is fundamentally non- 
perturbative in the strength of the disorder and enables 
us to study the RG flow of entire disorder distribu- 
tions, whereas the perturbative 1-loop momentum-shell 
RG based on the disorder averaged replica theory is re- 
stricted to the mean and the variance of the random mass 
distribution [22]. Both approaches, however, show that 
the relative variance diverges in the BG. This is an impor- 
tant result as it demonstrates that contrary to the gen- 
eral belief, the onset of Griffith instabilities is captured 
in perturbative RG. The comparison of the correlation- 
length exponents obtained by the two methods suggests 
that corrections beyond 1-loop order are small. 

The real-space RG approach presented here has a wide 
range of future applications. It can be used to study the 
effects of spatial correlations in the disorder and entails 
the search for self-similar disorder characterized by scale 
invariant distribution functions. Finally, the method is 
not restricted to disorder distributions, but can be used 
to study other inhomogeneities such as the so-called wed- 
ding cake structures of alternating MI and SF regions 
found in optical-lattice system [2"9"j . 

Acknowledgment: FK benefited from stimulating 
discussions with A. V. Chubukov, I. F. Herbut, J. Keel- 
ing, and R. Moessner. AH thanks B. Brinkman for many 
useful discussions. This work was supported by EP- 
SRC under grant code EP/I 004831/1 and by NSF-DMR- 
1104909. 



5 



M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. 
Fisher, Phys. Rev. B 40, 546 (1989). 
M. Greiner, O. Mandel, T. Hansch, and I. Bloch, Nature 
419, 51 (2002). 

R. T. Scalettar, G. G. Batrouni, and G. T. Zimanyi, 
Phys. Rev. Lett. 66, 3144 (1991). 

W. Krauth, N. Trivedi, and D. Ceperley, Phys. Rev. Lett. 
67, 2307 (1991). 

R. V. Pai, R. Pandit, H. R. Krishnamurthy, and S. Ra- 
masesha, Phys. Rev. Lett. 76, 2937 (1996). 
J. Kisker and H. Rieger, Phys. Rev. B 55, R11981 (1997). 
P. Sen, N. Trivedi, and D. M. Ceperley, Phys. Rev. Lett. 
86, 4092 (2001). 

J.-W. Lee, M.-C. Cha, and D. Kim, Phys. Rev. Lett. 87, 
247006 (2001). 

N. Prokof'ev and B. Svistunov, Phys. Rev. Lett. 92, 
015703 (2004). 

K. G. Singh and D. S. Rokhsar, Phys. Rev. B 46, 3002 
(1992). 

J. Wu and P. Phillips, Phys. Rev. B 78, 014515 (2008). 
U. Bissbort and W. Hofstetter, Europhys. Lett. 86, 50007 
(2009). 

R. Mukhopadhyay and P. B. Weichman, Phys. Rev. Lett. 
76, 2977 (1996). 

J. K. Freericks and H. Monien, Phys. Rev. B 53, 2691 
(1996). 



[15 
[16 
[17 
[18 

[19; 

[20 

[21 

[22 

[23 

[24 

[25 
[26 

[27 

[28; 
[29 



B. V. Svistunov, Phys. Rev. B 54, 16131 (1996). 

I. F. Herbut, Phys. Rev. Lett. 79, 3502 (1997). 

I. F. Herbut, Phys. Rev. B 57, 13729 (1998). 

P. B. Weichman and R. Mukhopadhyay, Phys. Rev. B 

77, 214516 (2008). 

F. Kriiger, J. Wu, and P. Phillips, Phys. Rev. B 80, 
094526 (2009). 

S. Iyer, D. Pekker, and G. Rafael, Phys. Rev. B 85, 
094202 (2012). 

L. Pollet, N. V. Prokof'ev, B. V. Svistunov, and 
M. Troyer, Phys. Rev. Lett. 103, 140402 (2009). 
F. Kriiger, S. Hong, and P. Phillips, Phys. Rev. B 84, 
115118 (2011). 

A. Aharony and A. B. Harris, Phys. Rev. Lett. 77, 3700 
(1996). 

F. Pazmandi, R. T. Scalettar, and G. T. Zimanyi, Phys. 

Rev. Lett. 79, 5130 (1997). 

A. B. Harris, J. Phys. C. 7, 1671 (1974). 

J. T. Chayes, L. Chayes, D. S. Fisher, and T. Spencer, 

Phys. Rev. Lett. 57, 2999 (1986). 

S. Sachdev, Quantum phase transitions (Cambridge 
University Press, Cambridge New York, 2011), ISBN 
0521514681. 

K. Sengupta and N. Dupuis, Phys. Rev. A 71, 033629 
(2005). 

N. Gemelke, X. Zhang, C.-L. Hung, and C. Chin, Nature 
460, 995 (2009). 



